home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_i.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  8.6 KB  |  306 lines
  1. /* specfunc/bessel_i.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_pow_int.h>
  26. #include <gsl/gsl_sf_bessel.h>
  27.  
  28. #include "error.h"
  29.  
  30. #include "bessel.h"
  31.  
  32.  
  33. /* i_{l+1}/i_l
  34.  */
  35. static
  36. int
  37. bessel_il_CF1(const int l, const double x, const double threshold, double * ratio)
  38. {
  39.   const int kmax = 2000;
  40.   double tk   = 1.0;
  41.   double sum  = 1.0;
  42.   double rhok = 0.0;
  43.   int k;
  44.  
  45.   for(k=1; k<=kmax; k++) {
  46.     double ak = (x/(2.0*l+1.0+2.0*k)) * (x/(2.0*l+3.0+2.0*k));
  47.     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
  48.     tk  *= rhok;
  49.     sum += tk;
  50.     if(fabs(tk/sum) < threshold) break;
  51.   }
  52.  
  53.   *ratio = x/(2.0*l+3.0) * sum;
  54.  
  55.   if(k == kmax)
  56.     GSL_ERROR ("error", GSL_EMAXITER);
  57.   else
  58.     return GSL_SUCCESS;
  59. }
  60.  
  61.  
  62. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  63.  
  64. int gsl_sf_bessel_i0_scaled_e(const double x, gsl_sf_result * result)
  65. {
  66.   double ax = fabs(x);
  67.  
  68.   /* CHECK_POINTER(result) */
  69.  
  70.   if(ax < 0.2) {
  71.     const double eax = exp(-ax);
  72.     const double y = ax*ax;
  73.     const double c1 = 1.0/6.0;
  74.     const double c2 = 1.0/120.0;
  75.     const double c3 = 1.0/5040.0;
  76.     const double c4 = 1.0/362880.0;
  77.     const double c5 = 1.0/39916800.0;
  78.     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
  79.     result->val = eax * sum;
  80.     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  81.   }
  82.   else if(ax < -0.5*GSL_LOG_DBL_EPSILON) {
  83.     result->val = (1.0 - exp(-2.0*ax))/(2.0*ax);
  84.     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  85.   }
  86.   else {
  87.     result->val = 1.0/(2.0*ax);
  88.     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  89.   }
  90.   return GSL_SUCCESS;
  91. }
  92.  
  93.  
  94. int gsl_sf_bessel_i1_scaled_e(const double x, gsl_sf_result * result)
  95. {
  96.   double ax = fabs(x);
  97.  
  98.   /* CHECK_POINTER(result) */
  99.  
  100.   if(ax < 3.0*GSL_DBL_MIN) {
  101.     UNDERFLOW_ERROR(result);
  102.   }
  103.   else if(ax < 0.25) {
  104.     const double eax = exp(-ax);
  105.     const double y  = x*x;
  106.     const double c1 = 1.0/10.0;
  107.     const double c2 = 1.0/280.0;
  108.     const double c3 = 1.0/15120.0;
  109.     const double c4 = 1.0/1330560.0;
  110.     const double c5 = 1.0/172972800.0;
  111.     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
  112.     result->val = eax * x/3.0 * sum;
  113.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  114.     return GSL_SUCCESS;
  115.   }
  116.   else {
  117.     double ex = exp(-2.0*ax);
  118.     result->val = 0.5 * (ax*(1.0+ex) - (1.0-ex)) / (ax*ax);
  119.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  120.     if(x < 0.0) result->val = -result->val;
  121.     return GSL_SUCCESS;
  122.   }
  123. }
  124.  
  125.  
  126. int gsl_sf_bessel_i2_scaled_e(const double x, gsl_sf_result * result)
  127. {
  128.   double ax = fabs(x);
  129.  
  130.   /* CHECK_POINTER(result) */
  131.  
  132.   if(ax < 4.0*GSL_SQRT_DBL_MIN) {
  133.     UNDERFLOW_ERROR(result);
  134.   }
  135.   else if(ax < 0.25) {
  136.     const double y = x*x;
  137.     const double c1 = 1.0/14.0;
  138.     const double c2 = 1.0/504.0;
  139.     const double c3 = 1.0/33264.0;
  140.     const double c4 = 1.0/3459456.0;
  141.     const double c5 = 1.0/518918400.0;
  142.     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
  143.     const double pre = exp(-ax) * x*x/15.0;
  144.     result->val = pre * sum;
  145.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  146.     return GSL_SUCCESS;
  147.   }
  148.   else {
  149.     double ex = exp(-2.0*ax);
  150.     double x2 = x*x;
  151.     result->val = 0.5 * ((3.0+x2)*(1.0-ex) - 3.0*ax*(1.0+ex))/(ax*ax*ax);
  152.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  153.     return GSL_SUCCESS;
  154.   }
  155. }  
  156.  
  157.  
  158. int gsl_sf_bessel_il_scaled_e(const int l, double x, gsl_sf_result * result)
  159. {
  160.   double sgn = 1.0;
  161.   double ax  = fabs(x);
  162.  
  163.   if(x < 0.0) {
  164.     /* i_l(-x) = (-1)^l i_l(x) */
  165.     sgn = ( GSL_IS_ODD(l) ? -1.0 : 1.0 );
  166.     x = -x;
  167.   }
  168.  
  169.   if(l < 0) {
  170.     DOMAIN_ERROR(result);
  171.   }
  172.   else if(x == 0.0) {
  173.     result->val = 0.0;
  174.     result->err = 0.0;
  175.     return GSL_SUCCESS;
  176.   }
  177.   else if(l == 0) {
  178.     gsl_sf_result il;
  179.     int stat_il = gsl_sf_bessel_i0_scaled_e(x, &il);
  180.     result->val = sgn * il.val;
  181.     result->err = il.err;
  182.     return stat_il;
  183.   }
  184.   else if(l == 1) {
  185.     gsl_sf_result il;
  186.     int stat_il = gsl_sf_bessel_i1_scaled_e(x, &il);
  187.     result->val = sgn * il.val;
  188.     result->err = il.err;
  189.     return stat_il;
  190.   }
  191.   else if(l == 2) {
  192.     gsl_sf_result il;
  193.     int stat_il = gsl_sf_bessel_i2_scaled_e(x, &il);
  194.     result->val = sgn * il.val;
  195.     result->err = il.err;
  196.     return stat_il;
  197.   }
  198.   else if(x*x < 10.0*(l+1.5)/M_E) {
  199.     gsl_sf_result b;
  200.     int stat = gsl_sf_bessel_IJ_taylor_e(l+0.5, x, 1, 50, GSL_DBL_EPSILON, &b);
  201.     double pre   = exp(-ax) * sqrt((0.5*M_PI)/x);
  202.     result->val  = sgn * pre * b.val;
  203.     result->err  = pre * b.err;
  204.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  205.     return stat;
  206.   }
  207.   else if(l < 150) {
  208.     gsl_sf_result i0_scaled;
  209.     int stat_i0  = gsl_sf_bessel_i0_scaled_e(ax, &i0_scaled);
  210.     double rat;
  211.     int stat_CF1 = bessel_il_CF1(l, ax, GSL_DBL_EPSILON, &rat);
  212.     double iellp1 = rat * GSL_SQRT_DBL_MIN;
  213.     double iell      = GSL_SQRT_DBL_MIN;
  214.     double iellm1;
  215.     int ell;
  216.     for(ell = l; ell >= 1; ell--) {
  217.       iellm1 = iellp1 + (2*ell + 1)/x * iell;
  218.       iellp1 = iell;
  219.       iell   = iellm1;
  220.     }
  221.     result->val  = sgn * i0_scaled.val * (GSL_SQRT_DBL_MIN / iell);
  222.     result->err  = i0_scaled.err * (GSL_SQRT_DBL_MIN / iell);
  223.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  224.     return GSL_ERROR_SELECT_2(stat_i0, stat_CF1);
  225.   }
  226.   else if(GSL_MIN(0.29/(l*l+1.0), 0.5/(l*l+1.0+x*x)) < 0.5*GSL_ROOT3_DBL_EPSILON) {
  227.     int status = gsl_sf_bessel_Inu_scaled_asymp_unif_e(l + 0.5, x, result);
  228.     double pre = sqrt((0.5*M_PI)/x);
  229.     result->val *= sgn * pre;
  230.     result->err *= pre;
  231.     return status;
  232.   }
  233.   else {
  234.     /* recurse down from safe values */
  235.     double rt_term = sqrt((0.5*M_PI)/x);
  236.     const int LMAX = 2 + (int) (1.2 / GSL_ROOT6_DBL_EPSILON);
  237.     gsl_sf_result r_iellp1;
  238.     gsl_sf_result r_iell;
  239.     int stat_a1 = gsl_sf_bessel_Inu_scaled_asymp_unif_e(LMAX + 1 + 0.5, x, &r_iellp1);
  240.     int stat_a2 = gsl_sf_bessel_Inu_scaled_asymp_unif_e(LMAX     + 0.5, x, &r_iell);
  241.     double iellp1 = r_iellp1.val;
  242.     double iell   = r_iell.val;
  243.     double iellm1 = 0.0;
  244.     int ell;
  245.     iellp1 *= rt_term;
  246.     iell   *= rt_term;
  247.     for(ell = LMAX; ell >= l+1; ell--) {
  248.       iellm1 = iellp1 + (2*ell + 1)/x * iell;
  249.       iellp1 = iell;
  250.       iell   = iellm1;
  251.     }
  252.     result->val  = sgn * iellm1;
  253.     result->err  = fabs(result->val)*(GSL_DBL_EPSILON + fabs(r_iellp1.err/r_iellp1.val) + fabs(r_iell.err/r_iell.val));
  254.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  255.  
  256.     return GSL_ERROR_SELECT_2(stat_a1, stat_a2);
  257.   }
  258. }
  259.  
  260.  
  261. int gsl_sf_bessel_il_scaled_array(const int lmax, const double x, double * result_array)
  262. {
  263.   int ell;
  264.   gsl_sf_result r_iellp1;
  265.   gsl_sf_result r_iell;
  266.   int stat_0 = gsl_sf_bessel_il_scaled_e(lmax+1, x, &r_iellp1);
  267.   int stat_1 = gsl_sf_bessel_il_scaled_e(lmax,   x, &r_iell);
  268.   double iellp1 = r_iellp1.val;
  269.   double iell   = r_iell.val;
  270.   double iellm1;
  271.   result_array[lmax] = iell;
  272.   for(ell = lmax; ell >= 1; ell--) {
  273.     iellm1 = iellp1 + (2*ell + 1)/x * iell;
  274.     iellp1 = iell;
  275.     iell   = iellm1;
  276.     result_array[ell-1] = iellm1;
  277.   }
  278.   return GSL_ERROR_SELECT_2(stat_0, stat_1);
  279. }
  280.  
  281.  
  282. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  283.  
  284. #include "eval.h"
  285.  
  286. double gsl_sf_bessel_i0_scaled(const double x)
  287. {
  288.   EVAL_RESULT(gsl_sf_bessel_i0_scaled_e(x, &result));
  289. }
  290.  
  291. double gsl_sf_bessel_i1_scaled(const double x)
  292. {
  293.   EVAL_RESULT(gsl_sf_bessel_i1_scaled_e(x, &result));
  294. }
  295.  
  296. double gsl_sf_bessel_i2_scaled(const double x)
  297. {
  298.   EVAL_RESULT(gsl_sf_bessel_i2_scaled_e(x, &result));
  299. }
  300.  
  301. double gsl_sf_bessel_il_scaled(const int l, const double x)
  302. {
  303.   EVAL_RESULT(gsl_sf_bessel_il_scaled_e(l, x, &result));
  304. }
  305.  
  306.